/* M Orsi - Essex group - Soton Univ - 2009 
   Direct evaluation of Lipid Lateral Diffusion (hence the name latDiff.c), mainly from Rapaport p. 124.
   Sets of diffusion measurements are begun at fixed time intervals, so that several 
   sets of measurements based on different time origins will be in progress simultaneously
   because of the overlapped measurements.
   Memory is allocated in module startUp.c, function AllocArrays().
*/

// compile: pgcc -o latDiff latDiff.c -lm

/* program include files */
#include "dataStructs.h"
#include "definitions.h"

/* C-library include files */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/* external data references */
extern const VecRTwo *xyLipidCom;
extern const VecR region; 
extern const int nLipids, stepCount;
extern const real deltaT, timeNow;
extern LatDiffBuff *latDiffBuff; /* structure LatDiffBuf defined in module dataStructs.h */
extern real *latSqrDiffAv;
extern int upperLayerVelCount, lowerLayerVelCount, 
  countLatDiffAv; //  I SHOULD LOCALISE THIS VARIABLE HERE
extern const int  limitLatDiffAv, /* total #individual estimates used to produce the average value of D */ 
  nBuffLatDiff, /* #sets of data being collected at any time, each occupying a separate storage buffer */
  nValLatDiff, /* #measurements contributing to the set used to produce a single unaveraged estimate of D */
  stepLatDiff; /* measurement frequency */
extern VecRTwo velSumUpperLayer, velSumLowerLayer;

static void PrintLatDiff()
{
  real tVal;
  int j;

  FILE *fp;
  char latDiffFile[ 20 ];
  char time_ns[ 9 ]; // string containing timeNow, e.g. "100", in ns

  //  sprintf( time_ns, "%d", stepCount ); // 
  sprintf( time_ns, "%d", Nint( timeNow *  TIME_ns ) );
  latDiffFile[ 0 ] = '\0';
  strcat( latDiffFile, "latDiff-" );
  strcat( latDiffFile, time_ns );
  strcat( latDiffFile, ".dat" );
  fp = fopen( latDiffFile, "w" );


  for ( j = 0; j < nValLatDiff; j++ ) {
    tVal = j * stepLatDiff * deltaT;
    //fprintf( fp, "%17.4f %15.10f\n", tVal * TIME_ps, latSqrDiffAv[ j ] * Sqr( LENGTH_CM ) / TIME_s * 1e5); /* [ps], [10^(-5)cm^2/s] */
    //fprintf( fp, "%17.4f %15.10f\n", tVal * TIME_ps, latSqrDiffAv[ j ] * Sqr( LENGTH_CM ) / TIME_s * 1e8); /* [ps], [10^(-8)cm^2/s] */
    //    fprintf( fp, "%12.5f %12.3f\n", tVal * TIME_ns, latSqrDiffAv[ j ] * Sqr( LENGTH_CM ) / TIME_s * 1e8); /* [ns], [10^(-8)cm^2/s] */
    fprintf( fp, "%12.5f %12.3f\n", tVal * TIME_ns, latSqrDiffAv[ j ] * Sqr( LENGTH_nm ) / TIME_mus ); /* [ns], [nm^2/microsecond] */
    //  fflush( fp );
  }
  fclose( fp );
}

static void ZeroLatDiff()
{
  int j;

  countLatDiffAv = 0;
  for ( j = 0; j < nValLatDiff; j++ ) latSqrDiffAv[ j ] = 0.;  /* resetting the calculation */
}


/* InitLatDiff. Initialises the calculation. Called from SetUpJob(). 
   The values initially assigned to "count" determine the spacing between measurements: negative initial values delay
   the start of data collection for each set of measurements until the appropriate moment. */
void InitLatDiff()
{
  int nb;

  for ( nb = 0; nb < nBuffLatDiff; nb++ ) latDiffBuff[ nb ].count = - nb * nValLatDiff / nBuffLatDiff; 
  ZeroLatDiff();
  VZeroTwo( velSumUpperLayer );
  VZeroTwo( velSumLowerLayer );
  upperLayerVelCount = lowerLayerVelCount = 0;
}

static void AccumLatDiff()
{
  real fac;
  int j, nb;
  
  for ( nb = 0; nb < nBuffLatDiff; nb++ ) {
    if ( latDiffBuff[ nb ].count == nValLatDiff ) { /* if buffer "nb" is full */
      for ( j = 0; j < nValLatDiff; j++ ) { /* scan through the buffer elements */
	latSqrDiffAv[ j ] += latDiffBuff[ nb ].latSqrDisp[ j ]; /* accumulate LatDiff */
      }
      latDiffBuff[ nb ].count = 0;  /* reset buffer counter */
      ++countLatDiffAv;
      if ( countLatDiffAv == limitLatDiffAv ) { /* if the #individual estimates required has been reached */
        fac = 1. / ( 4. * nLipids * stepLatDiff * deltaT * limitLatDiffAv ); // 4 <- (2 * D), being D = 2 in this case (lipid diffusion in the xy plane) 
        for ( j = 1; j < nValLatDiff; j++ ) {
	  latSqrDiffAv[ j ] *= fac / j;
	}
	//	printf("*** PrintLatDiff()...\n");
	PrintLatDiff();
	//	printf("*** PrintLatDiff()... done\n");
	//	printf("*** ZeroLatDiff()...\n");
	ZeroLatDiff();
	//	printf("*** ZeroLatDiff()... done\n");
      }
    }
  }
  //  printf("*** Exiting AccumLatDiff()... \n");
}

/* Called by SingleStep() every "stepLatDiff" steps (after equilibration, obviously). */
void EvalLatDiff()
{
  VecRTwo latDisp; // lateral lipid displacement (current - initial) in the xy plane 
  int n, nb, ni;

  for ( nb = 0; nb < nBuffLatDiff; nb++ ) { /* scan through the individual buffers */
    if ( latDiffBuff[ nb ].count == 0 ) { /* if empty buffer */
      DO_LIPID { // set origin vectors
        latDiffBuff[ nb ].orgTwo[ n ] = xyLipidCom[ n ];
        latDiffBuff[ nb ].latDispTrueTwo[ n ] = xyLipidCom[ n ];
      }
    }
    if ( latDiffBuff[ nb ].count >= 0 ) {
      ni = latDiffBuff[ nb ].count;
      latDiffBuff[ nb ].latSqrDisp[ ni ] = 0.;
      DO_LIPID {
        VSubTwo( latDisp, latDiffBuff[ nb ].latDispTrueTwo[ n ], xyLipidCom[ n ] );
        VDivTwo( latDisp, latDisp, region );
        latDisp.x = Nint( latDisp.x );
        latDisp.y = Nint( latDisp.y );
        VMulTwo( latDisp, latDisp, region );
        VAddTwo( latDiffBuff[ nb ].latDispTrueTwo[ n ], xyLipidCom[ n ], latDisp );
        VSubTwo( latDisp, latDiffBuff[ nb ].latDispTrueTwo[ n ], latDiffBuff[ nb ].orgTwo[ n ] );
        latDiffBuff[ nb ].latSqrDisp[ ni ] += VLenSqTwo( latDisp );
      }
    }
    ++latDiffBuff[ nb ].count;  /* increase buffer counter */
  }
  AccumLatDiff();
}

void main()
{
  if ( fp = fopen( xySoluteUnwrapped.dat, "r" ) == NULL ) {
    printf( "File xySoluteUnwrapped.dat could not be opened. Exiting...\n" ); Exit(0);
  }
  
}
